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Abstract 

The dynamics of energy relaxation in thermalized one- and two- 
dimensional arrays with nonlinear interactions depend in detail on 
the interactions and, in some cases, on dimensionality. We describe 
and explain these differences for arrays of the Fermi-Pasta-Ulam type. 
In particular, we focus on the roles of harmonic contributions to the 
interactions and of breathers in the relaxation process. 
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1 Introduction 



The ability of extended systems to support the localization and transport of 
vibrational energy has been invoked in a number of physical settings including 
DNA molecules |1] , hydrocarbon structures , energy storage and transport 
in proteins [§, |], ||, the creation of vibrational intrinsic localized modes 
in anharmonic crystals using an optimally chosen sequence of femtosecond 
laser pulses || , photonic crystal waveguides |7j , and targeted energy transfer 
between donors and acceptors in biomolecules ||. It has become increasingly 
clear that thermal fluctuations may strongly affect (sometimes leading to 
degradation but at times actually helping) the process of energy localization 
and energy mobility @, |, g 0, 0, [H], [TJ, [H], |T|, |T|, [T7|. It is thus clearly 



important to investigate the nature and dynamics of thermal fluctuations in 
nonlinear arrays. The understanding of the spatial and temporal evolution of 
the thermal relaxation landscape will in turn lead to a better understanding 
of other chemical and physical processes that may be occurring in the relaxing 
landscape. 

The spatial and temporal evolution of an influx of energy into or an ef- 
flux of energy out of a nonlinear array, and the dynamical pathways that 
characterize energy relaxation in such an array, depend on a large number of 
factors (for recent reviews of a huge literature see JL8|, [19|]). The nature of 
the interactions and where the nonlinearities reside (in the local potentials 
or in the interactions), the boundary conditions (free, fixed, or periodic), 
the size of the system, its dimensionality, the way in which energy is de- 
posited in the system (initial conditions), etc., can all influence the evolution 
profoundly, and no general formalism that encompasses all these variations 
has yet been developed. It is thus necessary in this study (as in most oth- 
ers) to circumscribe the range of our inquiry. We have been particularly 
interested in discrete extended systems in which localized energy can also 



be mobile [fQ , I3| , [14| , and so our studies have focused on arrays with hard 



nonlinear interactions and no local potentials, specifically on Fermi-Pasta- 
Ulam (FPU) lattices. Some movement of localized energy may also occur in 
arrays with soft local potentials (see e.g. [l|]), but it is much more difficult to 
achieve. We thus concentrate on FPU arrays. 



Specifically in this paper we study energy relaxation in one-dimensional [20 



and two-dimensional FPU arrays with quartic potentials. The Hamiltonian 
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in one dimension is 



j=l i=l i=l 

where iV is the number of sites; k and k' are the harmonic and anharmonic 
force constants, respectively. In two dimensions for aniVxiV lattice 

•^i— l,j) "I - {.'I'i.j - \ ) ] 

Xi-ij) 4 + (x i>3 - - Xij_i) 4 ] . (2) 

To study energy relaxation we initially thermalize the system at tempera- 
ture T (see below), then connect the boundary sites (two end sites for a Id 
system, 4(N — 1) edge sites for 2d arrays) to a zero-temperature reservoir 
via damping terms, and observe the thermal relaxation of the array toward 
zero temperature [ffO, EH], We find that relaxation occurs through 
an interesting cascade of decay times that is sensitively dependent on the 
precise form of the interactions. This leads us to focus on two issues in the 
relaxation process: 1) the explicit role of the harmonic terms vs. anharmonic 
terms in the FPU Hamiltonian, and 2) the effects of array dimensionality. 

To thermalize the system to a given temperature T we augment the equa- 
tions of motion resulting from Eqs. ([!]) with the Langevin prescription con- 
necting each site to a heat bath. In one dimension 

d 

£i = - — [V{xi - Xi-i) + V(x i+1 - x^] - 7 ii + r]i(t). (3) 

OXi 

Here V{xi — Xj) is the FPU potential, 70 is the dissipation parameter, and 
the rji(t) are mutually uncorrelated zero-centered Gaussian ^-correlated fluc- 
tuations that satisfy the fluctuation-dissipation relation at temperature T: 

(rk(t)) = 0, (Vi{t)Vj{t')) = 2 lo k B TS lJ 5(t - t'). (4) 

The brackets here and below denote an ensemble average, and k B is Boltz- 
mann's constant. The generalization to two dimensions is immediate. We 
implement free-end boundary conditions, that is, xq = X\ and xn = xn+i in 



N ■ 2 , N 

i,j=l i,j=l 

V N 
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one dimension and x j = xij, x^,j = x^+xj, %i,o — x i,i, an d x i,N — x i,N+i in 
two dimensions. For the integrations of the equations of motion we use the 
fourth order Runge-Kutta method. 

The thermalization process involves the spontaneous emergence of anhar- 
monic (including localized) modes. In order to understand the dynamics of 
such modes, we begin again with a thermalized array and explicitly inject a 
breather-like excitation of energy much higher than the thermal energy. We 
then observe how the entire system, thermalized array plus injected excita- 
tion, evolves and relaxes toward zero temperature. 

In Sec. we itemize the measures used to display the relaxation process. 
In Sec. [| we present our relaxation results, obtained primarily from numerical 
simulations. Section | shows what happens to an explicitly injected local 
excitation as the arrays relax to zero temperature. Section || provides a brief 
summary and synthesis of the outcomes. 

2 Measures of Thermal Relaxation 

Once our system is thermalized to temperature T we disconnect it from the 
thermal bath (i.e., we remove the i]i(t) and 7o£; terms from Eq. (H) or the 
equivalent terms from the two-dimensional set of equations) , and we connect 
the edge sites to a cold T = reservoir through the addition of dissipative 
terms — 72 to the equations of motion for these sites. We then continue the 
integration using the thermalized positions and displacements as the initial 
conditions. An ensemble is constructed by repeating this experiment for 
different thermalized initial conditions. There are of course a variety of ways 
to display the outcome, and in our work we have chosen three. One, the 
most "global" measure, is to follow the decay of the total array energy E(t) 
as a function of time. The second is to follow the spectrum of the system 
as a function of time. This gives a frequency-by-frequency picture of the 
relaxation process and therefore a more complete description. The third is 
to show a spatial energy landscape as a function of time. 

The total energy E(t) of the array is simply the Hamiltonian function 
evaluated with the displacements and velocities obtained from the equations 
of motion. For a one- dimensional harmonic array this function has been 
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calculated analytically by Piazza et al. [|TU| for small 7 and large N: 
E{t) 1 



- rdqe-^=e-^I (t/r ) 

7T Jo 



E(0) J 



e -t/r for t <^T , 

(2£)-" 2 fort>>To . < 5 > 

In the first line Jo is the modified zero-order Bessel function, and tq = N/2 r j. 
The decay time r(q) for phonons of wavevector q is 

1 2 (q 



fees 2 r- . (6) 



r(q) r \2 

The second line in Eq. (|5|) gives the short-time and long-time behaviors. The 
former is a simple exponential decay associated with the lowest frequency 
phonon mode since it has the shortest decay time. The power law relaxation 
arises from the cascade of different decay times of the different phonon modes. 
For finite N the decay becomes exponential again when only the modes near 
the band-edge of the spectrum still survive. Note that the decay is not expo- 
nential throughout (a common misapprehension for harmonic systems), al- 
though the initial exponential behavior does last longer the larger the system. 
Each phonon mode decays separately (exponentially) and independently of 
all the others. This translates to an independent and separate decay of each 
frequency portion of the spectrum (see below). The calculation of the total 
energy is more cumbersome but basically similar for two-dimensional arrays. 
One finds that the decay time for phonons of wavevector q = (q x ,q y ) is 



1 



2 fQx\ . 2 fly 
cos — + cos — 
V2/ V2 



which upon integration over wavevectors leads to 



(7) 



E(t) , \e- 2t/T0 fort<r , 

^ = ^('/-) = { (¥) -' fori>>T0 . <*> 

The behavior of E(t)/E(0) for Id and 2d anharmonic arrays is expected 
to be different than Eqs. fl5|) and (§), respectively. These behaviors will be 
presented in the next section. 
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Our second measure of thermal relaxation focuses on the decay pathways 
of the different spectral regions as the array cools down. We concentrate 
on the time evolution of the Fourier transform of the relative displacement 
correlation function. Relative displacements provide a particularly sensitive 
measure of how adjacent masses are moving relative to one another. In a 
thermalized array the spectrum of interest is defined as 

POO 

S(u) = 2 dr C(r) cos cut, (9) 
Jo 

where in one dimension 

1 N 

C(r) = J2 <W* + r)<W*)> (10) 

and SiXi(t) is the relative displacement 

5iXi{t) = Xi(t) - Xi-xit). (11) 

In two dimensions 

j N N 

C(r) = _ E ( SiXi ^ + T )^ X M) 

\ I j=2 j=X 



+ MiV _ n E E fo**^ + r )<^M w> (12) 

V ' 8=1 j=2 

where 

SiXij(t) = X id (t) - Xi-x,j(t), SjX id (t) = X id (t) - SCij-].^). (13) 



The thermal equilibrium spectrum for harmonic arrays in one and two 
dimensions can be calculated analytically. In one dimension with periodic 
boundary conditions (for sufficiently long chains the boundary conditions do 
not affect the equilibrium spectrum) 

_ 4-f k B T ^ 1 - cos(27rg/iV) 
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where 



n, 2 ( 5 ) = -! ±v /(!) 2 -«si„*(f). (15) 



In two dimensions 

JV-l 

2 - COS[2TTp/l\) - COS{2TT 

r o / \. on r o / \. 

■ uo 



= 4 l0 k B T ^ 2-cos(27rp/iV)-cos(27rg/iV) 



^ 2 [r 2 (p,g)+^][r 2 (p,g) + 

where now 



7o , //7o 



2 V V 2 
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ri l2 (p, ?) = -£ ± A/ ( £ ) -4* sin 2 ( ^) + sin 2 ( ^ ) . (17) 



NJ \N/\ 



For anharmonic chains these spectra must be obtained numerically. 

To monitor the decay of the spectrum when the thermalized arrays are 
connected to a cold reservoir we introduce the time-dependent spectra 

PTmax 

S(u,t) = 2 dr C(r,t) cos lot (18) 

Jo 

where r max = 2n/uo m i n and LO m i n is chosen for a desired frequency resolution; 
the choice L0 min = 0.0982, corresponding to r max = 64, turns out to be 
numerically convenient. The time-dependent correlation function is actually 
an average over the time interval t — to to t, where we have chosen t — 100 
(short enough for the correlation function not to change appreciably but long 
enough for statistical purposes) and is defined as follows in one dimension: 

C ^ *) = Jn^T) S a* J dT ' {6lXl{t " T ' )6tXl{t ~ T '~ T)) ' (19) 

where At = to — r max . The generalization to the two-dimensional case is 
obvious. 



3 Thermal Relaxation 

We start by presenting two sets of figures, each associated with one of the 
measures for relaxation mentioned in the previous section. Each set presents 
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Figure 1: Temporal relaxation of the fractional energy for different arrays 
with potential parameters and initial temperatures indicated in the panels. 
First panel: one dimension (N = 50). Second panel: two dimensions (20 x 20 
lattices). In all cases 7 = 0.1. The fractional energy for the harmonic arrays 
is independent of temperature. 

both Id and 2d results. Along with these results we present some auxiliary 
figures that help in the interpretation of the outcomes. 

The total energy of a relaxing array decays in time, and the questions 
of interest are how exactly the energy decreases with time for arrays with 
different interactions and in different dimensions. The answers are illustrated 
in Fig. III. Accompanying these decay curves are the more detailed spectral 
decay curves shown in Fig. |2|. 

The decay of the total energy ratio for harmonic systems is independent 
of temperature, which is verified numerically and therefore leads to a single 
curve in Fig. [TJ in each dimension for the given parameters k, N, and 7. The 
initial exponential decay of the energy in both one and two dimensions, and 
the N and 7 dependences of the decay rate, have been verified numerically. 
The long-time decay as an inverse power law E(t)/E(0) ~ t~ a for both cases 
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has also been verified. It is interesting to note that Eq. (|5]) can be rewritten 
as an integral that makes explicit the cascade of relaxation times giving rise 
to the inverse power law behavior: 



-t/r 



t/T0 Ut/r ) = dT ; - , (20) 



E (°) * Jr /2 t^/t-to/2' 

The lower limit, due to a nonzero shortest decay time, leads to the initial 
exponential decay of the energy, but it is the long r~ 3 ^ 2 tail of slow relaxation 
times that leads to the power law decay. We return to this point below. 

The associated spectra for the harmonic arrays are shown in each of the 
first frames of the two panels in Fig. 0. The evolution confirms that low fre- 
quencies decay more rapidly in the harmonic chain - the spectrum is absorbed 
by the cold reservoir from the bottom up, and by the latest times shown, only 
the longer-lived band-edge modes remain in the system. Since each spectral 
component is associated with an independent phonon, the spectral decrease 
occurs "vertically" , that is, each spectral component decays directly into the 
reservoir on its characteristic time scale; this is shown schematically for the 
one-dimensional systems in Fig. [|, where the downward arrows represent ab- 
sorption by the reservoir and their relative length schematizes the absorption 
rate. 

For anharmonic arrays the relaxation behavior depends strongly on the 
presence or absence of a harmonic component in the potential, and, in some 
respects, on dimensionality. First we consider the purely anharmonic arrays. 
The dimensionality plays a particularly important role in this case. 

In one dimension the upper panel of Fig. [I] shows an essentially purely 
exponential decay (verified separately), which is characteristic of a single pre- 
dominant decay channel. The decay is more rapid at higher temperatures. 
Note that there are no phonons in this purely quartic system, so that single 
frequencies are not associated with normal modes of the system. Conversely, 
exact solutions of the FPU chain such as solitons and intrinsic localized 
modes may involve many frequencies. The second frame in the upper panel 
of Fig. H shows that the higher frequencies decay first, exactly opposite to the 
harmonic chain. We find that the relaxation pathway is for the high frequency 
portions of the spectrum to degrade rapidly into lower frequency excitations, 
as schematically indicated by the sloped arrows in Fig. |3]. The lowest fre- 
quencies decay into the reservoir and define the exponential decay rate seen 
in Fig. ffl. In more detail, the high frequency components of the spectrum are 
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mainly due to highly mobile localized modes that degrade into lower energy 
(less mobile) excitations as they move and collide with one another. The 
lowest frequency excitations are in turn absorbed into the cold reservoir but 
continue to be replenished through the degradation process. It is important 
to note, however, that among the low frequency excitations are some that 
persist for a very long time, certainly beyond the times of our simulations. 
Their decay is surely slower than exponential, perhaps a stretched exponen- 
tial. These, the only remaining spectral components at time t = 2000, are 
"labeled" by short downward arrows in the relaxation schematic and include 
rather stable breather and/or soliton modes that move very slowly and are 
localized away from the boundaries. Also, some direct relaxation of all fre- 
quency components into the reservoir occurs as well (shown schematically by 
the short arrows at high frequencies in Fig. ||), but this direct relaxation is 
slower than the energy degradation pathway. For example, when a highly 
mobile localized excitation reaches a boundary, it typically remains at the 
boundary for about one period of oscillation (which is short for a highly en- 
ergetic excitation), during which it loses a small portion of its energy to the 
reservoir. The remaining excitation is reflected back into the chain, where it 
will continue to lose energy through other collision events and/or re-arrival 
at the boundaries. The role of high-frequency mobile modes and of low- 
frequency slowly moving or stationary modes in this picture will be tested in 
more detail in the next section, where we explicitly inject a high-frequency 
localized mode into the array and observe the relaxation dynamics. We do 
note here that our picture is consistent with known facts about localized 
states. In particular, it is known that higher- frequency and/or higher ampli- 
tude localized modes can move at higher velocities [|| [1^, [TfJ]. It is also known 
that while in motion such modes lose energy through collisions with other 
excitations. Figure [I| shows a faster decay at higher temperatures, which 
is consistent with our observations elsewhere that the speed of an injected 
pulse (and therefore, we conjecture, the speed of a moving localized mode) 
in these arrays increases with temperature fl3|| . 

The relaxation dynamics of the purely anharmonic array in two dimen- 
sions differs from the one- dimensional case in a number of ways. First, we 
note that the decay of the energy in Fig. |I] is (except for very early times) 
slower in the hard array than in the harmonic one over the times of obser- 
vation. Second, except for an initial short time interval, the decay is found 
not to be exponential, indicating that there is no longer a single predominant 
decay channel as there was in one dimension. In two dimensions the relax- 
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ation pathway again includes degradation of higher frequency excitations to 
lower frequencies, as can be seen in the spectral rendition in Fig. ^|. However, 
whereas in one dimension the degradation process is faster than the decay 
of low frequency excitations into the reservoir (and hence this latter decay is 
the rate-limiting step that defines the total energy decay process), here the 
degradation process is slower, leading to spectral bottlenecks and competing 
time scales. We find that increasing the array size leads to slower degra- 
dation of the high frequency components and to more pronounced spectral 
bottlenecks in the mid-frequency range. Our physical picture of the source 
of the competing time scales involves the observation that in two dimen- 
sions, localized excitations are not nearly as mobile as in one dimension. In 
one dimension energy degradation occurs as a consequence of high mobility 
and the resultant inevitable frequent scattering. The reduced mobility in 
two dimensions was noted in earlier work on pulse propagation O] and will 
be supported and clarified in the next section, where we explicitly inject a 
high-amplitude breather into our system. At long times in the 2d system 
excitations of all energies eventually reach the boundaries. These excitations 
typically lose some of their energy into the cold reservoir and the remainder 
is reflected back into the array, where it either degrades into lower energy 
components or reaches a boundary again with the attendant energy loss. 
Note also that with increasing temperature the total system energy decays 
more rapidly, which is consistent with our assertion that mobility (low as 
it may be) in the purely hard arrays increases with temperature because of 
the participation of more mobile higher frequency components in the initial 
equilibrium mix of excitations. 

The result is an energy decay that is of stretched exponential form, 



as a function of time for four purely anharmonic 2d arrays of different sizes, 
damping coefficients, and temperatures. When the decay is a stretched ex- 
ponential this rendition gives a flat line at the value (3{t) = a. We find that 
the value of a is independent of temperature, as seen in the figure (we have 





as is evident in the first panel of Fig. |], where we plot 
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tested this assertion for various lattice sizes in the range T = 0.1 — 1.0), and 
that it increases with 7, as also seen in the figure. For the 20 x 20 lattices 
we find that a is around 0.33 for 7 = 0.1 and a ~ 0.43 for 7 = 1. The 
temperature independence is explained by the fact that temperature does 
not dominate the mobility of the residual excitations nor does it determine 
their rate of energy loss once they reach a boundary. On the other hand, it 
is reasonable that increasing 7 leads to a faster long-time relaxation process 
(larger a) because more energy is lost to the reservoir upon each collision. 
We find that a also increases with N, which is somewhat of a puzzle. We 
find, for instance, that in a 50 x 50 array with 7 = 0.1 (and for any temper- 
ature) a ~ 0.52. The issue of the size dependence of relaxational processes 
following a stretched exponential behavior is a difficult problem that has 



only recently been addressed in a different context |£3] . In summary, for the 
purely anharmonic 2d arrays 

a = a(iV,7). (23) 

The mixed arrays, i.e., those with interactions that have both quadratic 
and quartic potential contributions, are of course the ubiquitous FPU systems 
since it is difficult to envision a "real" physical system that has no quadratic 
potential terms (one might also say this about cubic potential terms that 
have not been included here pi]). The thermal relaxation of these arrays 



proceeds similarly in one and two dimensions. At early and intermediate 
times the mixed arrays relax very similarly to the harmonic arrays (i.e., ex- 
ponential decay followed by power law decay), albeit somewhat modified and 
speeded up by the presence of high-frequency mobile excitations in addition 
to the low-frequency phononic excitations. The rapid decay of both low and 
high frequency excitations is evident in the third frames of both panels in 
Fig. |2|. This similarity in one dimension was noted by Piazza et al. fTUfl . 
After some time, however, the mixed chain relaxation behavior changes to a 
stretched exponential. This occurs when the low frequency modes have es- 
sentially all decayed. The higher frequency spectral components that persist 
are localized long-lived excitations (note that high-frequency phonon modes 



are unstable against breather formation in these systems [12]). Unlike the 
purely anharmonic array, here there is no low-frequency residue to perturb 
the high-frequency localized excitations and so they survive relatively un- 
perturbed and immobile for a long time. The persistence of high-frequency 
spectral components is seen in the third frames of both panels in Fig. [| and 
the schematic representation of the progression in one dimension is shown 
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in Fig. [5| The slow leakage of breather energy into low energy modes that 
continue to dissipate into the cold reservoir is responsible for the eventual 
stretched exponential relaxation of the system. 

These behaviors can be seen clearly in the second panel of Fig. [| which 
shows j3(t) for various Id and 2d mixed arrays. The initial exponential be- 
havior (which corresponds to a flat curve at (3(t) = 1), occurs over too short 
a time scale to be clearly discernible. There then follows a power law decay 
regime which eventually turns to a stretched exponential. In the power law 
regime, if the energy decays as E(t)/E(0) ~ (t/r ) _a then it follows that 

(3{t) = ^ln , which is independent of a and depends only on the ratio 

iV/7 via tq. The two 2d curves and the Id curve with the same value of iV/7 
are all seen to decay with the same slope. As all the curves settle into their 
asymptotic behavior, we see that a depends neither on N nor on 7: the two 
Id arrays are of different sizes but asymptote to the same a (approximately 
0.21), as do the two 2d arrays of the same size and initial temperature but 
with different damping coefficients [a ~ 0.03). On the other hand, the 2d 
arrays that have different initial temperatures asymptote to different values 
of a (approximately 0.25 for T = 0.1 and around 0.03 for T = 0.5). The 
observed behavior is explained by the fact that the slow leakage of energy 
out of the long-lived localized excitations is rate limiting. A higher initial 
temperature leads to more energetic, more stable breathers with slower leak- 
age, hence explaining why a decreases with increasing temperature. Neither 
the size of the system nor the damping coefficient are important in this limit, 
since the slowest process is the leakage. The low-energy outcome of that 
leakage is absorbed quickly by the cold reservoir. Note that this description 
is consistent with that provided for relaxation of lattices with local anhar- 



monic potentials [22]. In summary, for Id and 2d arrays with quadratic plus 



quartic interaction potentials 

a = a(T). (24) 

It is interesting to stress that, like an inverse power law decay, a stretched 
exponential can indeed be obtained from a distribution or hierarchical pro- 
gression of decay times. For example, 
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which should be compared with Eq. (POD- More generally, if the distribution 
of relaxation times varies as e _ ^ T//Tl - )M , then the decay will go as a stretched 
exponential with o ~ 1/(1 + \i) at long times (the stretched exponential 
with a = 1/2 is the only one whose associated distribution is expressible 
in extremely simple analytic form, but distributions associated with other 
fractional exponents are also known analytically or numerically [0). The 
distributions leading to inverse power decay and to a stretched exponential 
decay are both broad, but the inverse power law of course arises explicitly 
from very long tails not present in the stretched exponential. In other words, 
the relaxation of the last energy residues of a very large harmonic lattice 
take longer than those of a very large anharmonic lattice. There is an impor- 
tant difference, however: in the harmonic lattice the persistent excitations 
are distributed over the entire lattice, while in the anharmonic lattice they 
are localized. Also, in a finite lattice eventually the harmonic decay is again 
exponential while that of the anharmonic system remains a stretched expo- 
nential. 

Spatial energy landscapes rendered in gray scale provide a pictorial rep- 
resentation of the thermalization progressions. We follow widespread con- 
vention and define the local energy in one dimension as 

E * = Y + l ~ Xi -^ + v ^ Xl+1 ~ (26) 

with obvious generalization in two dimensions. Figure [| shows the temporal 
evolution of the local energy landscape for each of the three chains. Par- 
ticularly dramatic is the spontaneous occurrence of an essentially stationary 
breather in the mixed chain. It is this sort of breather that leads to the 
extremely slow relaxation of the mixed chain energy. Spatial landscapes for 
two-dimensional systems are shown in six time frames in Figs. ^ and |7| for 
the pure anharmonic and the mixed anharmonic lattices, respectively. Note 
that in the purely hard array the localized high-energy regions move around 
and relax within the time scale of the progression. In the mixed array, on 
the other hand, the "hot spots" persist and essentially do not move. 

4 Relaxation With an Injected Breather 

Our portrayal of the thermal relaxation process can be further bolstered 
by initially injecting a high-energy localized excitation in the center of each 
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thermalized array and observing the behavior of this excitation during the 
relaxation process. Some of the dimensionality differences are thereby clar- 
ified. The injected excitation is chosen so as to be close to a known exact 
breather solution of the anharmonic arrays. 

In a one-dimensional array with an interaction potential V(xi — Xi-i) = 
(xi — Xi_i) n as n — ► oo and exact odd-parity breather is one of amplitude 
A at a site and —A/2 at each of the two immediately adjacent sites. An 
exact even-parity breather is one with amplitude A at one site and —A at 
an immediately adjacent site. These are not exact solutions when n is not 
infinite and/or when there are quadratic contributions to the potential, but 



they are close to exact, even for the FPU chain |26[ P7| . In two dimensions 
the odd-parity solution with amplitude A at one site and amplitudes —A/ 4 
at the four nearest neighboring sites is also nearly exact, but there is no 
equivalent to the even-parity breather. We insert an odd-parity excitation 
at the center of our array and in each case choose arrays sufficiently large 
(300 sites in one dimension, 30 x 30 in two dimensions) so that the excitation 
either decays or stops moving before ever reaching the boundaries. In har- 
monic arrays the fate of such an injected excitation is completely predictable 
and uninteresting: it spreads quickly over the entire array and thus loses 
its localized character. The associated Fourier decomposition into phonon 
modes dictates the relaxation behavior. 

To follow the excitation in the one- dimensional anharmonic arrays we 
calculate the mean squared displacement 

<^(*)> = ((W(*) - y) 2 ) (27) 

as a measure of the position of the excitation (its dispersion in the anharmonic 
chains is very small Here N/2 is the initial point of highest energy in 

the chain and i ma x{t) is the point of maximum energy at time t. Similarly, 
in two dimensions we define 

(r 2 {t)) = (ii X!max {t) - j) 2 ^ + /(z tf|Tnaa! (t) - y) 2 ^ (28) 

where (N/2, N/2) is the initial point of highest energy and (i x ,max(t),iy,max(t)) 
is the point of maximum energy at time t. 

In a purely anharmonic chain in one dimension, a short time after it is 
created the excitation begins to move essentially ballistically in one direction 
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or the other with equal probability. The motion continues for a period of 
random duration, until the excitation stops moving for a random period of 
time. Then it moves again in either direction. Whatever its initial parity, 
while subsequently stationary the excitation has even parity (the more stable 
of the two configurations). Any perturbation (usually scattering of slow low- 
frequency excitations) that disturbs this parity sets the breather in motion, 
and while it moves it alternates between even and odd parity. The excitation 
only loses energy while in motion, through collisions with persistent low- 
frequency excitations. 

A typical mean squared displacement for the one- dimensional purely an- 
harmonic array is shown in Fig. §. The mean square displacement follows 
the superdiffusive law (x 2 (t)) ~ t a with a = 1.5 over the entire lifetime of the 
excitation. This particular exponent is recovered for the purely quartic chain 
under all conditions that we have tested, that is, independently of force con- 
stant, excitation amplitude, and temperature. Parameter variations affect 
only the prefactor, which reflects the breather velocity. Indeed, it does not 
matter when in the course of the relaxation process the localized excitation 
is introduced: its mean squared displacement grows with the same exponent 
1.5 until the excitation is extinguished into the background. This confirms 



the role of the persistent low-frequency excitations. We have argued [20 



that this particular superdiffusive exponent can be understood in terms of 
scattering events that occur with a time distribution v{t) ~ t~ 5 / 2 [28 . 



In a purely anharmonic two-dimensional array the motion is rather dif- 
ferent (which is consistent with the spectral differences in one and two di- 
mensions). A typical mean squared displacement is also shown in Fig. || 
The law is now sub diffusive, (r 2 (t)) ~ t a with a = 0.89 over the lifetime 
of the excitation. Again, this exponent is insensitive to force constant, ex- 
citation amplitude, and temperature changes. It does reflect the fact that 
the excitation moves much less in two dimensions than in one (although it 
does of course move and eventually recedes into the background). In ear- 
lier work [13| we noted that a pulse in a one dimensional purely hard array 
tends to move more rapidly but remains more tightly concentrated than a 
pulse in, say, a harmonic or soft array. We also noted that in two (or more) 
dimensions these two tendencies are in some sense contradictory since the 
only way that a symmetric excitation can move is by breaking its symmetry 
and/or dispersing. The sort of perturbation that would set a breather in 
motion requires an asymmetry that is more difficult to achieve in two dimen- 
sions than in one. If the distribution of collision times of energetic breathers 
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with other excitations that can set it in motion has sufficiently long quiescent 
periods (long-tailed waiting time probability distribution function) then the 
motion of the excitations is typically subdiffusive [p^ . Note that this does 
not preclude collisions that lead to energy loss by the breather even if it is 
not set in motion. 

The situation in mixed anharmonic arrays is similar in one and two di- 
mensions, see Fig. |8|. The injected excitations at first move with the same 
characteristic exponents a as in the corresponding purely anharmonic sys- 
tems, but as the harmonic interactions sweep the background thermal energy 
out of the system, the localized excitations stop moving. The mean squared 
displacement in each dimension then becomes independent of time (earlier 
in two dimensions than in one). 



5 Summary 

Energy relaxation in one- and two-dimensional nonlinear arrays with quartic 
interparticle interactions (Fermi- Past a- Ulam or FPU arrays) proceeds along 
energetic pathways completely different from those of harmonic systems and 
is quite sensitive to the presence or absence of quadratic contributions to the 
interactions. Relaxation in a purely harmonic array involves the sequential 
decay of independent phonon modes starting with those of lowest frequency 
and moving upward across the spectrum. The decay of energy in these arrays 
is exponential at short and long times, but follows an inverse power law at 
intermediate times. Throughout the decay process, the energy is distributed 
uniformly over the entire array. A localized excitation injected in the lattice 
simply decays according to the distribution of characteristic relaxation times 
of its phonon components. 

FPU arrays with quadratic and quartic interactions contain phonon-like 
modes as well as high-energy nonlinear, and to varying degrees localized, ex- 
citations (provided the initial temperature is sufficiently high to excite these). 
The relaxation process involves the decay of phonons in the same spectral 
order as in the harmonic arrays and also energy losses through collisions of 
mobile high-energy localized excitations as they collide with lower-frequency 
ones. Eventually the harmonic interactions succeed in "sweeping" the system 
clean of low energy excitations and the remaining localized modes are quasis- 
tationary breather solutions that persist for a very long time. The decay to 
this quasistationary state is a stretched exponential with an exponent that 
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depends on temperature but not on system size or damping coefficient. An 
explicitly injected high-energy localized breather follows this behavior, with 
a mean squared displacement that at first grows with time but then becomes 
independent of time when the system has been swept clean of background 
excitations. 

The greatest differences between one- and two-dimensional FPU arrays 
occur in the purely quartic arrays. In all cases the relaxation begins from the 
high-frequency end of the spectrum (opposite to the harmonic case) and in- 
volves not a direct absorption by the cold reservoir but rather a degradation 
of high-frequency excitations to lower-frequency ones. In one dimension this 
degradation process is considerably faster than in two dimensions. The low- 
est frequency modes in one dimension are absorbed by the cold reservoir but 
are quickly replenished by the degradation process. The energy relaxation is 
essentially exponential, with a time constant determined by the decay of the 
lowest frequency components. In two dimensions the degradation process is 
slower and there are frequency bottlenecks so that the decay of the lowest 
frequency excitations into the cold reservoir no longer constitutes the rate 
limiting process. Instead, the degradation and decay contribute to a result- 
ing stretched exponential energy decay with an exponent that depends on 
system size and damping coefficient but is independent of temperature. In 
both one and two dimensions there remains a thermal residue of localized 
low-frequency excitations that continue to perturb and degrade higher fre- 
quency ones. The array is never "swept clean" of low-energy excitations as is 
the mixed array, and therefore no persistent breathers occur in this system. 
In order to confirm this behavior we have followed the dynamics of an in- 
jected high-frequency localized excitation in these arrays. In one dimension 
this localized excitation remains localized but is very mobile throughout its 
lifetime, being characterized by a super- diffusive mean squared displacement. 
In two dimensions the excitation also remains localized and is much less mo- 
bile (but nevertheless, always mobile until it disappears), being characterized 
by sub-diffusive motion. 
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Figure 2: Time evolution of spectra for various relaxing arrays. The upper 
three-frame panel is for one-dimensional (50 sites) arrays, the lower three- 
frame panel for two-dimensional (20 x 20 sites) arrays. First frames: har- 
monic interactions (k = 0.5). Second frames: purely anharmonic interactions 
(k' = 0.5). Third frames: mixed interactions (k — k' — 0.5). The time pro- 
gression is as indicated. The t — spectrum (solid curves) in each case is 
the equilibrium spectrum at T = 0.5, the initial temperature. In all cases 
7 = 0.1. The thin vertical lines indicate the frequencies oo = \^4k = V2 (Id 
panels) and u = 2 (2d panels). 
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Figure 3: Schematic representation of the spectral relaxation channels in 
one dimension. The one-dimensional spectra of Fig. ^ for times t = and 
t = 2000 are shown again here, and the arrows depict the pathways of differ- 
ent spectral components. Downward arrows indicate absorption by the cold 
reservoir, while angled arrows denote degradation from one spectral region 
to another. The relative lengths of the arrows depict the associated rates. 
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Figure 4: Plot of (3(t) as a function of time. A flat line below (3{t) = 1 
indicates stretched exponential behavior. First panel: purely anharmonic 
2d arrays of different sizes, damping coefficients, and temperatures. Second 
panel: various Id and 2d mixed arrays. 
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Figure 5: Local energy landscapes of one-dimensional 30-site arrays initially 
thermalized at T = 0.5. Time advances along the y-axis until t = 1000. A 
gray scale is used to represent the local energy, with darker shading corre- 
sponding to more energetic regions. First panel: harmonic chain, k = 0.5 
and k! = 0. Second panel: purely anharmonic chain, k — and k! = 0.5. 
Third panel: mixed chain, k — k' — 0.5. 



24 



Figure 6: Local energy landscapes for a 20 x 20 purely anharmonic lattice 
(k' = 0.5) initially thermalized at T = 0.5 and with 7 = 0.1. From first to 
last frames: * = 0, 400, 800, 1200, 1600, and 2000. 
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Figure 7: Local energy landscapes for a 20 x 20 mixed anharmonic lattice 
(k = k' = 0.5) initially thermalized at T = 0.5 and with 7 = 0.1. From first 
to last frames: t = 0, 400, 800, 1200, 1600, and 2000. 
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Figure 8: Mean squared displacement of a localized mode in various one 
dimensional (300 sites) and two dimensional (30 x 30) FPU arrays. The 
slopes of the two short straight lines are 1.5 and 0.89. In all cases A = 2, 
T = 0.1, and 7 = 0.1. 
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